home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dsytrd.f < prev    next >
Text File  |  1997-06-25  |  9KB  |  280 lines
  1.       SUBROUTINE DSYTRD( UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO )
  2. *
  3. *  -- LAPACK routine (version 2.0) --
  4. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  5. *     Courant Institute, Argonne National Lab, and Rice University
  6. *     September 30, 1994
  7. *
  8. *     .. Scalar Arguments ..
  9.       CHARACTER          UPLO
  10.       INTEGER            INFO, LDA, LWORK, N
  11. *     ..
  12. *     .. Array Arguments ..
  13.       DOUBLE PRECISION   A( LDA, * ), D( * ), E( * ), TAU( * ),
  14.      $                   WORK( * )
  15. *     ..
  16. *
  17. *  Purpose
  18. *  =======
  19. *
  20. *  DSYTRD reduces a real symmetric matrix A to real symmetric
  21. *  tridiagonal form T by an orthogonal similarity transformation:
  22. *  Q**T * A * Q = T.
  23. *
  24. *  Arguments
  25. *  =========
  26. *
  27. *  UPLO    (input) CHARACTER*1
  28. *          = 'U':  Upper triangle of A is stored;
  29. *          = 'L':  Lower triangle of A is stored.
  30. *
  31. *  N       (input) INTEGER
  32. *          The order of the matrix A.  N >= 0.
  33. *
  34. *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  35. *          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
  36. *          N-by-N upper triangular part of A contains the upper
  37. *          triangular part of the matrix A, and the strictly lower
  38. *          triangular part of A is not referenced.  If UPLO = 'L', the
  39. *          leading N-by-N lower triangular part of A contains the lower
  40. *          triangular part of the matrix A, and the strictly upper
  41. *          triangular part of A is not referenced.
  42. *          On exit, if UPLO = 'U', the diagonal and first superdiagonal
  43. *          of A are overwritten by the corresponding elements of the
  44. *          tridiagonal matrix T, and the elements above the first
  45. *          superdiagonal, with the array TAU, represent the orthogonal
  46. *          matrix Q as a product of elementary reflectors; if UPLO
  47. *          = 'L', the diagonal and first subdiagonal of A are over-
  48. *          written by the corresponding elements of the tridiagonal
  49. *          matrix T, and the elements below the first subdiagonal, with
  50. *          the array TAU, represent the orthogonal matrix Q as a product
  51. *          of elementary reflectors. See Further Details.
  52. *
  53. *  LDA     (input) INTEGER
  54. *          The leading dimension of the array A.  LDA >= max(1,N).
  55. *
  56. *  D       (output) DOUBLE PRECISION array, dimension (N)
  57. *          The diagonal elements of the tridiagonal matrix T:
  58. *          D(i) = A(i,i).
  59. *
  60. *  E       (output) DOUBLE PRECISION array, dimension (N-1)
  61. *          The off-diagonal elements of the tridiagonal matrix T:
  62. *          E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
  63. *
  64. *  TAU     (output) DOUBLE PRECISION array, dimension (N-1)
  65. *          The scalar factors of the elementary reflectors (see Further
  66. *          Details).
  67. *
  68. *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
  69. *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  70. *
  71. *  LWORK   (input) INTEGER
  72. *          The dimension of the array WORK.  LWORK >= 1.
  73. *          For optimum performance LWORK >= N*NB, where NB is the
  74. *          optimal blocksize.
  75. *
  76. *  INFO    (output) INTEGER
  77. *          = 0:  successful exit
  78. *          < 0:  if INFO = -i, the i-th argument had an illegal value
  79. *
  80. *  Further Details
  81. *  ===============
  82. *
  83. *  If UPLO = 'U', the matrix Q is represented as a product of elementary
  84. *  reflectors
  85. *
  86. *     Q = H(n-1) . . . H(2) H(1).
  87. *
  88. *  Each H(i) has the form
  89. *
  90. *     H(i) = I - tau * v * v'
  91. *
  92. *  where tau is a real scalar, and v is a real vector with
  93. *  v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in
  94. *  A(1:i-1,i+1), and tau in TAU(i).
  95. *
  96. *  If UPLO = 'L', the matrix Q is represented as a product of elementary
  97. *  reflectors
  98. *
  99. *     Q = H(1) H(2) . . . H(n-1).
  100. *
  101. *  Each H(i) has the form
  102. *
  103. *     H(i) = I - tau * v * v'
  104. *
  105. *  where tau is a real scalar, and v is a real vector with
  106. *  v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i),
  107. *  and tau in TAU(i).
  108. *
  109. *  The contents of A on exit are illustrated by the following examples
  110. *  with n = 5:
  111. *
  112. *  if UPLO = 'U':                       if UPLO = 'L':
  113. *
  114. *    (  d   e   v2  v3  v4 )              (  d                  )
  115. *    (      d   e   v3  v4 )              (  e   d              )
  116. *    (          d   e   v4 )              (  v1  e   d          )
  117. *    (              d   e  )              (  v1  v2  e   d      )
  118. *    (                  d  )              (  v1  v2  v3  e   d  )
  119. *
  120. *  where d and e denote diagonal and off-diagonal elements of T, and vi
  121. *  denotes an element of the vector defining H(i).
  122. *
  123. *  =====================================================================
  124. *
  125. *     .. Parameters ..
  126.       DOUBLE PRECISION   ONE
  127.       PARAMETER          ( ONE = 1.0D0 )
  128. *     ..
  129. *     .. Local Scalars ..
  130.       LOGICAL            UPPER
  131.       INTEGER            I, IINFO, IWS, J, KK, LDWORK, NB, NBMIN, NX
  132. *     ..
  133. *     .. External Subroutines ..
  134.       EXTERNAL           DLATRD, DSYR2K, DSYTD2, XERBLA
  135. *     ..
  136. *     .. Intrinsic Functions ..
  137.       INTRINSIC          MAX
  138. *     ..
  139. *     .. External Functions ..
  140.       LOGICAL            LSAME
  141.       INTEGER            ILAENV
  142.       EXTERNAL           LSAME, ILAENV
  143. *     ..
  144. *     .. Executable Statements ..
  145. *
  146. *     Test the input parameters
  147. *
  148.       INFO = 0
  149.       UPPER = LSAME( UPLO, 'U' )
  150.       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  151.          INFO = -1
  152.       ELSE IF( N.LT.0 ) THEN
  153.          INFO = -2
  154.       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  155.          INFO = -4
  156.       ELSE IF( LWORK.LT.1 ) THEN
  157.          INFO = -9
  158.       END IF
  159.       IF( INFO.NE.0 ) THEN
  160.          CALL XERBLA( 'DSYTRD', -INFO )
  161.          RETURN
  162.       END IF
  163. *
  164. *     Quick return if possible
  165. *
  166.       IF( N.EQ.0 ) THEN
  167.          WORK( 1 ) = 1
  168.          RETURN
  169.       END IF
  170. *
  171. *     Determine the block size.
  172. *
  173.       NB = ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 )
  174.       NX = N
  175.       IWS = 1
  176.       IF( NB.GT.1 .AND. NB.LT.N ) THEN
  177. *
  178. *        Determine when to cross over from blocked to unblocked code
  179. *        (last block is always handled by unblocked code).
  180. *
  181.          NX = MAX( NB, ILAENV( 3, 'DSYTRD', UPLO, N, -1, -1, -1 ) )
  182.          IF( NX.LT.N ) THEN
  183. *
  184. *           Determine if workspace is large enough for blocked code.
  185. *
  186.             LDWORK = N
  187.             IWS = LDWORK*NB
  188.             IF( LWORK.LT.IWS ) THEN
  189. *
  190. *              Not enough workspace to use optimal NB:  determine the
  191. *              minimum value of NB, and reduce NB or force use of
  192. *              unblocked code by setting NX = N.
  193. *
  194.                NB = MAX( LWORK / LDWORK, 1 )
  195.                NBMIN = ILAENV( 2, 'DSYTRD', UPLO, N, -1, -1, -1 )
  196.                IF( NB.LT.NBMIN )
  197.      $            NX = N
  198.             END IF
  199.          ELSE
  200.             NX = N
  201.          END IF
  202.       ELSE
  203.          NB = 1
  204.       END IF
  205. *
  206.       IF( UPPER ) THEN
  207. *
  208. *        Reduce the upper triangle of A.
  209. *        Columns 1:kk are handled by the unblocked method.
  210. *
  211.          KK = N - ( ( N-NX+NB-1 ) / NB )*NB
  212.          DO 20 I = N - NB + 1, KK + 1, -NB
  213. *
  214. *           Reduce columns i:i+nb-1 to tridiagonal form and form the
  215. *           matrix W which is needed to update the unreduced part of
  216. *           the matrix
  217. *
  218.             CALL DLATRD( UPLO, I+NB-1, NB, A, LDA, E, TAU, WORK,
  219.      $                   LDWORK )
  220. *
  221. *           Update the unreduced submatrix A(1:i-1,1:i-1), using an
  222. *           update of the form:  A := A - V*W' - W*V'
  223. *
  224.             CALL DSYR2K( UPLO, 'No transpose', I-1, NB, -ONE, A( 1, I ),
  225.      $                   LDA, WORK, LDWORK, ONE, A, LDA )
  226. *
  227. *           Copy superdiagonal elements back into A, and diagonal
  228. *           elements into D
  229. *
  230.             DO 10 J = I, I + NB - 1
  231.                A( J-1, J ) = E( J-1 )
  232.                D( J ) = A( J, J )
  233.    10       CONTINUE
  234.    20    CONTINUE
  235. *
  236. *        Use unblocked code to reduce the last or only block
  237. *
  238.          CALL DSYTD2( UPLO, KK, A, LDA, D, E, TAU, IINFO )
  239.       ELSE
  240. *
  241. *        Reduce the lower triangle of A
  242. *
  243.          DO 40 I = 1, N - NX, NB
  244. *
  245. *           Reduce columns i:i+nb-1 to tridiagonal form and form the
  246. *           matrix W which is needed to update the unreduced part of
  247. *           the matrix
  248. *
  249.             CALL DLATRD( UPLO, N-I+1, NB, A( I, I ), LDA, E( I ),
  250.      $                   TAU( I ), WORK, LDWORK )
  251. *
  252. *           Update the unreduced submatrix A(i+ib:n,i+ib:n), using
  253. *           an update of the form:  A := A - V*W' - W*V'
  254. *
  255.             CALL DSYR2K( UPLO, 'No transpose', N-I-NB+1, NB, -ONE,
  256.      $                   A( I+NB, I ), LDA, WORK( NB+1 ), LDWORK, ONE,
  257.      $                   A( I+NB, I+NB ), LDA )
  258. *
  259. *           Copy subdiagonal elements back into A, and diagonal
  260. *           elements into D
  261. *
  262.             DO 30 J = I, I + NB - 1
  263.                A( J+1, J ) = E( J )
  264.                D( J ) = A( J, J )
  265.    30       CONTINUE
  266.    40    CONTINUE
  267. *
  268. *        Use unblocked code to reduce the last or only block
  269. *
  270.          CALL DSYTD2( UPLO, N-I+1, A( I, I ), LDA, D( I ), E( I ),
  271.      $                TAU( I ), IINFO )
  272.       END IF
  273. *
  274.       WORK( 1 ) = IWS
  275.       RETURN
  276. *
  277. *     End of DSYTRD
  278. *
  279.       END
  280.